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Abstract 

We study a calcium aluminosilicate glass of composition (Si02)o.67-(Al203)o.i2- (CaO)o.2i by 
means of molecular-dynamics (MD) simulations, using a potential made of two-body and three- 
body interactions. In order to prepare small samples that can subsequently be studied by first- 
principles, the finite size effects on the liquid dynamics and on the glass structural properties are 
investigated. We find that finite size effects affect the Si-O-Si and Si-O-Al angular distributions, 
the first peaks of the Si-0, Al-0 and Ca-0 pair correlation functions, the Ca coordination and the 
oxygen atoms environment in the smallest system (100 atoms). We give evidence that these finite 
size effects can be directly attributed to the use of three-body interactions. 

PACS numbers: 75.40.Mg,61.43.Fs,31.15.Gg 
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I. INTRODUCTION 



Materials used in nuclear wastes confinement are made of glasses with complex com- 
positions, and are based mainly on an alumino-boro-silicate network. These glasses have 
been studied for decades with the aim of testing their durability and their chemical resis- 
tance to water leaching, the process believed to be the limiting factor for their long time 
durabilit y^i^i^i^i^i^i'^'i^i^ . When such glasses are altered by water leaching, the layers formed 
on the surface are hydrated silica-like porous gels enriched with aluminum and calcium. 
Although it is often argued that the protective property of this layer stops the alteration 
processes induced by the leaching, no real justification has yet been giveniSJiiSiHiii. 

In calcium aluminosilicate glasses (CAS), silicon and aluminum are both network form- 
ers and have a tetrahedral coordination with oxygeusia. The tetrahedra are connected by 
shared oxygen atoms which are called bridging oxygens (BO) and they form an amorphous 
network in which Ca atoms are dispersed. Depending on the Al203/CaO ratio, Ca atoms 
play the role of modifiers if they create non-bridging oxygens (NBO) by breaking T-O-T 
linkages (T=Si/Al) and/or they play the role of charge balancing, by neutralizing the AIO4 
entities. In principle, the number of NBO atoms in a CAS glass is entirely determined by the 
Al203/CaO ratio and this number is directly related to the viscosity of the glass forming 
liquid. However it has been shown that this simple picture is not fully correct and that 
an excess number of NBO atoms exists even for an Al203/CaO ratio equal to one, which 
corresponds - in theory - to a perfect polymerized network. This excess number of NBO 
atoms could explain the viscosity anomaly measured in these systemsi^ and could also be 
related to the existence of other structural units such as oxygen triclusters and/or highly 
coordinated Si and Al atoms. The calcium environment is also not very well defined in that 
it is quite difficult to attribute a specific coordination to this species. Nevertheless, in pure 
CAS systems, it was observed by X-ray and neutron spectroscopies that the distribution of 
oxygen neighbors around Ca atoms is independent of the compositio n^^i^^ . Another feature 
observed in aluminosilicate glasses is the Al/Al avoidance principle: As the formation of 
Al-O-Al linkages is energetically less favorable than Si-O-Al linkages, the occurrence of the 
former should be very low in CAS struct ure o^^i^° . 

Due to the experimental difficulties to study thin layers, an important first step to obtain 
a microscopic understanding of these materials can be achieved by the use of molecular- 
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dynamics (MD) simulations. The aim of the present work is therefore to obtain a reasonable 
description of a CAS glass which composition is closed to the composition of the protective 
porous layer. However, in the case of such complex ternary systems, it is extremely difficult 
to find classical potentials that are able to give simultaneously a reliable description of the 
ionic and covalent bond lengths, of the bond angles and of the dynamics of the system. 
Alternatively, one can obtain a good description of such structural characteristics and of 
the dynamics of these complex systems from ab initio calculations. But in this case, the 
price to pay is the very limited size of the system (100 to 200 atoms) and the very short 
accessible length of the MD trajectories, thus making difficult the generation of disordered 
systems with good statistics on the one hand, and the equilibration of the samples at low 
temperatures on the other hand. In the aim of testing the classical potentials, we intend to 
use an approach in which the structure of a glass generated by classical MD is relaxed - and 
therefore tested - by ab initio MD. This approach has already been successfully used in the 
case of a Si02 glasst^i and of a Na20-4Si02 glass^^. Since ab initio simulations are limited 
to a rather small number of atoms, we first carried out a systematic study of the finite size 
effects on the aluminosilicate system of interest using a classical potential. 

In this paper, we present the study of a calcium aluminosilicate glass by MD simulations in 
which we analyzed the dependence of the structural and dynamical properties on the system 
size. In the first section, technical details on the MD simulations are given. Subsequently 
the dynamics of the liquids are presented in section two and the third section is devoted to 
the results obtained on the glass structure. Finally we discuss and summarize the results. 

II. METHODOLOGY 

The chosen composition of the calcium aluminosilicate glass samples for this study is 67% 
Si02 - 12% AI2O3 - 21% CaO in mole percent. 

To describe the interatomic interactions, we used Born-Mayer-Huggins (BMH) and 
Stillinger- Weber (SW)^^*2ii2^ potentials with the modified parameters of Delaye^^ which were 
obtained from an earlier study of complex multi-oxide glasses. The potentials include two- 
body and three-body terms. In Eq. the two-body interaction is given, where Vij denotes 
the distance between atoms i and j, and Aij, Btj, Cij and Dij are adjustable parameters. 
The first term of this pair potential corresponds to the Coulomb term, the second one to 
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the core repulsion, and the third and fourth ones are referred to as multipolar terms: 

F.,(r.,) = ^ + A., exp(-5., r.,) - (^) ' + (^) ' • (1) 

The Coulomb interactions were calculated using formal charges equal to +4 for Si, -2 for 
O, +3 for Al and +2 for Ca, and the standard Ewald summation method with alpha L=7.2 
was used where L is the size of the cubic box. Bond orientations in the CAS system can be 
described by adding a three-body term involving the 0-Si-O, 0-Al-O and Si-O-Si triplets. 
The three-body interaction is given by the product of a two-body interaction for each pair 
of the considered triplets by an angular strained force. 




-cos(^5J)^ (2) 



where r^j and rik denote the distances between atoms i and j (or atoms i and A;), Ojik is 
the angle between atoms j, i and k, and Ajjfc, 7ij, 7^^, and r^^, are adjustable parameters. 
Tables HI and ITTl summarize the parameters of the two-body and the three-body interactions. 

Newton's equations were integrated using a Velocity Verlet scheme with a time step of 
1.8 fs accepting a final deviation of +7 K on the total energy for MD runs of 1 million steps. 

Liquid samples of 100 and 200 atoms, generated using a random starting configuration, 
were equilibrated at different temperatures (from ^ 3000 K to ~ 7000 K). The lowest 
temperature for which the silicon self-diffusion was at least 7 A after 1.8 ns of MD trajectory 
for both system sizes, was taken as the liquid equilibration temperature: it corresponds to 
roughly 4200 K. Then, we generated and equilibrated independent liquid samples of different 
sizes at this temperature during 1.8 ns: 8 samples of 100 atoms, 4 samples of 200 atoms, 
2 samples of 400 atoms, 2 samples of 800 atoms and 1 sample of 1600 atoms. The sizes 
of the corresponding boxes were respectively 11.32 A, 14.26 A, 17.97 A, 22.64 A and 28.52 
A, corresponding to a density of 2.42 g.cm~^, which is the equilibrium density obtained 
with this potential at this temperature. The experimental density for glasses of similar 
compositions at 300 K lies between 2.55 and 2.60 g/cw?^^. 

III. DYNAMICS OF THE LIQUID 

The mean square displacements (MSD), {r'^if)) = ((rj(0) — riit))"^), of each atomic species 
were calculated at different temperatures for the systems of 100 and 200 atoms. In order 
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to obtain a good description of the dynamical properties, it appeared that an average over 
a larger number of samples was needed. Therefore these dynamical properties were studied 
only for the smallest system sizes of 100 and 200 atoms and using averages over 30 indepen- 
dent samples. Figure Q represents the MSD of the various types of atoms in the 200-atoms 
liquid function of time in a logarithmic plot. 

At short times, all curves exhibit a straight line with a slope of 2 which indicates that 
all atoms move ballistically. Then, after a transient regime, the atoms show a diffusive 
behavior which is characterized by straight lines with slopes of 1. From the mean squared 
displacements at long times, the self-diffusion constants D were extracted using Einstein's 
law: 

B=itam. (3) 

The resulting D of the different species for systems of 200 atoms as a function of the inverse 
temperature are depicted in Fig. |21 The most diffusive element is Ca, followed by O, Al 
and Si. It is interesting to note that the Al atoms follow the oxygen self-diffusion for all 
temperatures, despite the fact that, as Si, it is a network former. This behavior resembles the 
one found in the Al203-Si02 system^^. Only at the lowest temperature the Al self-diffusion 
approaches the one of the silicon atoms. The activation energies of the self- diffusion of the 
Si, O, Al and Ca atoms were extracted from a linear fit of the curves in Fig. [21 (for 1/T 
greater than 2.3 10^'^ K^^) and are equal to 6.4 eV, 5.9 eV, 6.5 eV and 5.3 eV, respectively. 

The evolution of the diffusion constants as a function of temperature as well as the 
values of the activation energies do not show any differences between the 100-atoms and 
200-atoms liquids, suggesting that no size effects exist - at least up to 200 atoms - for the 
dynamical quantities. This is in apparent contrast with the results of previous studies 
carried out on supercooled liquids^i^ and in silica^ which shows strong finite size effects 
in the dynamics of the smallest samples. However, in order to study more deeply the finite 
size effects on the dynamics of the CAS liquids, more statistics would have been needed for 
systems of larger sizes as well. Since our main purpose is to study the size effects on the 
glass structural properties, we considered that this is beyond the scope of the present work. 
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IV. STRUCTURAL ANALYSIS 



After cooling down all liquids from 4200 K to 300 K using a linear quench with a rate of 
10^^ K.s and after a MD simulation of 150 ps at 300 K, we performed 1000 additional 
steps for the structural analysis. 



Radial distribution functions 



The radial distribution functions (RDF) were evaluated by means of Eq.(@]) where Na 
and Nj3 denote the number of atoms for species a and f3 and V is the volume^: 



" ie{a}je{f3} 



Figure El shows the radial distribution functions for Si-0, Al-0 and Ca-0 pairs for all system 
sizes, with the 100-atoms system as dashed-lines. In all systems apart from the 100-atoms 
ones, the positions of the maxima of the Si-0 and Al-0 RDFs, which correspond to the 
most probable bond distances, are respectively equal to 1.60 A and 1.76 A. These results 
are in good agreement with the X-ray diffraction measurements by Petkov et alM^ on CAS 
glasses of close compositions which give values of 1.63 A and 1.60 A for Si-0 and 1.77 A and 
1.75 A for Al-0. Because of the complexity of the Ca environment, more statistics would 
have been needed in order to refine the location of the maximum of the Ca-0 RDF first 
peak. The latter can nevertheless be assigned approximately to 2.50 A. The lower panels 
in Fig. El present zooms on the Si-0 (a), Al-0 (b), and Ca-0 (c) RDF first peaks. In these 
zooms, one can clearly notice that the systems of 100 atoms exhibit a different RDF than 
the ones of larger sizes since the Si-0 and Al-0 first peaks are narrower and higher. This 
indicates that the local order is more pronounced in the 100-atoms samples. The zoom on 
the Ca-0 RDF clearly shows two peaks in the maximum. This double peak is not due to 
poor statistics but appears to be due to correlations that exist in all samples of 100 atoms. 

In Fig. 12 we compare the Si-Ca, Al-Ca, Al-Al and Ca-Ca RDFs for the systems of 100- 
(dashed lines) and 1600-atoms (bold lines). Due to the small number of Al and Ca atoms, 
the RDF curves are quite noisy for both system sizes. The determination of interatomic 
distances is consequently quite difficult, particularly for the Ca-Ca distance. For systems of 
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100-atoms, we can remark that the RDFs show several peaks that are absent in the larger 
system. Nevertheless, the RDFs of the 100-atoms system follow the global shape of the 
RDFs of the system with 1600 atoms, Fig. EJd. The same behavior is observed for the Si-Si, 
Si-Al and 0-0 RDFs^s. 

The first peak of the Si-Ca and Al-Ca RDFs of the 1600-atoms system presents a double 
shoulder. Fig. and Fig. which corresponds to two kinds of Si-Ca distances at 3.37 
A and 3.57 A. A more pronounced double shoulder is also present on the first peak of the 
Al-Ca RDFs at 3.37 A and 3.57 A . These double shoulders are due to the presence of non- 
bridging atoms (NBO) in the Si04 or AIO4 groups. As mentioned in the introduction, the 
non-bridging oxygen (NBO) atoms are defined as oxygen atoms that have only one network 
former (Si or Al) as nearest neighbor. Here, these neighbors are counted using cutoff radii 
defined by the location of the minima after the first peaks in the Si-0 and Al-0 RDFs. 

A detailed analysis of the Si-Ca RDF is presented in Fig. Et^, b and c. Figure ]^ shows 
the first peak of the Si-0 RDF for the different environments Q^^ of the Si atoms. Since a 
Qx atom is defined as a network former atom surrounded by x bridging oxygen atoms, a Qx 
tetra-coordinated network former is thus surrounded by (4-x) non-bridging oxygen atoms. 
Figure clearly shows that the first peak of the Si-0 RDF shifts to lower distances when 
the number of NBO is increased on the Si04 groups (i.e. x is decreased). The Si-0 bond 
length shifts from 1.60Atol.52A. Fi gureEb presents the Ca-0 RDFs by differentiating BO 
and NBO atoms. A shorter Ca-0 bond length is observed when the oxygen is non-bridging. 
By differentiating Q^; Si atoms in the Si-Ca RDFs, Fig. Eb, it can be clearly concluded that 
the first shoulder, observed in the first peak of the Si-Ca RDF, is due to NBO atoms and 
the second to BO atoms, since the separation of the two peaks becomes more and more 
pronounced as the Si atoms are more and more surrounded by NBO atoms. 

From the analysis of the 1600-atoms system, we summarize all interatomic distances of 
the CAS structure in Tab. IIIH the error bars being of the order of 5- 10"'^ A in the case 
of Si-0 and Al-0 and between 10~^ A and 5-10~^ A for the other distances. In the case 
of Si-Ca and Al-Ca, we took the positions of the first peak as a crude estimation of the 
distances, the RDF first peak being too noisy to allow a clear determination of the position 
of the maximum. Some of the values given in Tab. IIIII can be compared with experimental 
distances obtained either on glasses or on crystalline analogous. Interatomic distances do 
not change dramatically with composition and one can compare directly the Si-0 and the 
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Al-0 distances with the experimental values presented in Tab. IIIII 

In summary, we find clear finite size effects for the 100-atoms system which affect the 
first peaks of the Si-0, Al-0 and Ca-0 RDFs. These finite size effects suggest that the glass 
models become more locally organized for the smallest system size. For all the other system 
sizes, the RDFs give information on the local structure which is overall in agreement with 
experimental data except for the Ca-0 bond length that is found to be slightly too large. 

B. Angular distributions and linkages 

In order to calculate angular distributions, we have defined cutoff radii using the location 
of the minima after the first peaks in the Si-0 and Al-0 RDFs, Fig. EK- Figure IHl represents 
the Al-O-Al (a), Si-O-Si (b), Si-O-Al (c), 0-Al-O (d) and 0-Si-O (e) angular distributions 
for the systems of 100 and 1600 atoms. All other system sizes show similar distributions as 
the ones of the 1600-atoms system. Like the RDFs, the angular distributions of the systems 
of 100 atoms follow the general shape of the 1600-atoms ones but they present many small 
peaks that are due to the more ordered character of the systems. Particularly in the case of 
the Si-O-Si angular distribution, we have checked that these peaks are not due to statistics 
but correspond indeed to structural correlations that are present at the same position in all 
samples of 100 atoms. 

In the 1600-atoms case, the 0-Al-O and 0-Si-O angular distributions. Fig. IHli and e, show 
a quite symmetric shape centered at 107° and 108°, respectively, which is in agreement with 
the values measured in amorphous silicates^. The Si-O-Si and Si-O-Al angular distributions. 
Fig. inb and c, have maxima at 157° and 148° respectively, and exhibit an asymmetric shape 
with a tail toward small angles. The Si-O-Si angles are larger than the Si-O-Al ones which 
is consistent with the results of static ab initio calculations by Xiao and Lasaga^^ in which 
they showed that the equilibrium Si-O-Si angle of a H6Si20 cluster is larger than the Si-O-Al 
angle of a HgAlSiO cluster. The Al-O-Al angular distribution. Fig. IHt, spreads from 90 to 
180° and does not present a distinctive shape with a clear maximum like in the Si-O-Si and 
Si-O-Al distributions. Because of the low number of Al-O-Al linkages in the system, a better 
statistics would be needed in order to refine the Al-O-Al distribution. That a few Al-O-Al 
linkages are present in the glass samples is consistent with the fact that the Al/Al avoidance 
principle^^*^ is not necessarily fulfilled in amorphous aluminosilicates. 
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In order to estimate to what extent the Al/Al avoidance principle is violated in our 
glass samples, we computed the number of T-O-T' linkage concentrations (T, T' = Si or 
Al) that would be obtained if the Si and Al atoms were randomly distributed among all 
the possible network former sites in the 1600-atoms sample^^. We then compared these 
concentrations with those found in the simulated glasses. If the Si"'"'' and Al"*"^ cations 
were randomly distributed among the network former sites of the CAS glass, the number of 
Al-O-Al linkages would be given by : 

where A^aiai is the number of Al-O-Al linkages, A'^o is the number of bridging oxygens, 
A'sio and A^aio are the number of Si-0 and Al-0 bonds, respectively. The probability to 
find an Al-O-Al linkage corresponds to the probability to pick out two Al-0 bonds among 
the different possibilities for the Si-0 and Al-0 bonds. In the system of 1600 atoms, the 
expectation value for the number of Al-O-Al linkages obtained from Eq. © is 62 whereas 
only 39 Al-O-Al linkages are found in the simulated sample. This difference indicates that 
the Al/Al avoidance exists for this kind of glass but is not completely fulfilled. 

Fig. [7| shows the mean number of Si-O-Si, Si-O-Al, Al-O-Al linkages as well as the total 
number of linkages, divided by the total number of oxygen atoms, as a function of the 
system size. The variation of the number of linkages does not exceed 1 % and remains 
inside the error bars, thus indicating that within the statistical errors no finite size effects 
can be detected on this quantity and that size effects only appear on the associated angular 
distributions. 



C. Coordination numbers 



In order to determine the Si, Al and Ca coordination numbers, we used the same cutoff 
radii as the ones used for the angular distributions. We found that the local order of the Si 
and Al atoms is very well defined in that the atoms are always fourfold coordinated whatever 
the size of the system is. However, as the minimum of the Ca-0 RDF is not defined very 
well, we calculated the Ca coordination distributions for ten different cutoff radii, chosen so 
that they sample regularly the region of the minimum (see inset of Fig. |HJ. Only five of these 
distributions, as well as their average, are shown in Fig. |Hlfor the 1600-atoms system. The 
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coordination distribution of the Ca atoms is strongly dependent on the choice of the cutoff 
radii and therefore caution should be used while discussing the average Ca coordination 
number. 

The averaged Ca coordination distributions for the different system sizes are presented 
in Fig. El From this figure, we can see that the Ca atoms have multiple local environments, 
the number of oxygen surrounding the Ca atoms ranging from 4 to 12. The maxima of these 
distributions are located between 7 and 8 except for the 100-atoms systems. In that case, 
the distribution seems to be shifted toward smaller numbers, to be more asymmetric and to 
have a maximum located at 6. From experimental measurements, the coordination of the 
Ca atoms were estimated to be equal to 7 in WAXS^^ and to 5.2 in X-ray diffraction^i, and 
to 6 by Car-Parrinello simulations in a CAS liquid of 100 atoms^. Therefore it seems that 
our models slightly overestimate the Ca mean coordination and that this overestimation is 
compensated by finite size effects in the 100-atoms case. 

D. Bridging oxygens, non-bridging oxygens, and oxygen triclusters 

As defined in Sec. IIV Al the non-bridging oxygen (NBO) atoms are oxygen atoms that 
have only one network former (Si or Al) as nearest neighbor. These neighbors are counted 
using the cutoff radii used for the angular distributions as well as for the coordination 
numbers. 

The expected number of NBO atoms in the system can be estimated from simple stoi- 
chiometry arguments, if the system is supposed to be made of perfect tetrahedra with only 
2-fold oxygen atoms. In that case, the number of NBO atoms is simply given by: 

A^NBO = 2 X A^ca - A^Ai. (6) 

Following this argument, the percentage of NBO atoms in all systems should be equal to 
9.52 %, if the Ca^"*" ions fully play their charge balancing role. However, in the simulated 
samples, we always find an excess number of NBO atoms. Fig. ITUk shows the percentage 
of NBO as a function of the system size and this percentage is clearly located between 10.7 
% and 11.8 %. The excess number of NBO atoms could be due to the use of a too high 
quench rate, however recent ^^O NMR experiments on calcium and sodium aluminosilicate 
glassesi^i^ revealed the presence of excess NBO atoms as well. 
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One can observe a slight size dependence for this number indicating that, for larger 
systems, a larger number of NBO atoms is found. It is furthermore interesting to note that 
a large majority of NBO atoms are located on the Si tetrahedra and only a few percent 
of them on Al tetrahedra (between 1.7 and 8.0 % of the NBO atoms are located on Al 
tetrahedra). This result is in agreement with X-ray diffraction experiments by Petkov et 
al^ who observed that the NBO atoms are exclusively located on Si tetrahedra in CAS 
glasses. One can also estimate the number of NBO atoms that are located on a given Si 
tetrahedron, i.e. if the Si tetrahedron is a Q4 (no NBO), a Q3 (one NBO), a Q2 (two NBO) 
etc. This information is given in Tab. IIVI As it can be seen, there are no clear size effects 
on the relative number of Q^, and most of the Si tetrahedra with NBO atoms are in the Q3 
conformation. 

It has been proposed that the excess number of oxygen atoms could be compensated by 
oxygen triclusters, i.e. oxygen atoms which are bonded to three network former atomsi^jii. 
We therefore evaluated the number of oxygen triclusters in our samples and we found that 
the percentage varies between 1 % and 1.5 %, Fig. ITUh . Again, this percentage shows an 
increase with the system size as for the number of NBO atoms. Fig. ITUk . even if the error 
bars are quite large for the small system sizes. Finally the number of BO atoms (which are 
not triclusters) was also evaluated and the resulting percentage is depicted in Fig. ITUb . 

In order to determine if there exists any particular trend in the formation of tricluster 
types (OSis, 0Si2Al, 0SiAl2 and OAI3), we evaluated the percentage of these units in the 
case of a random distribution of the Si and Al atoms among the possible network former 
sites. Let A^^os be the number of oxygens triclusters. Then there are SxNq^ bonds between 
the oxygen triclusters and the three network former neighbors. Among these bonds, we can 
distinguish the number of Si-03 bonds, given by A^^si-os, and the number of Al-03 bonds, 
given by A^ai-03- We consider the entity formed by the oxygen tricluster and its three 
neighbors as an ensemble of three T-0 bonds (T = Si or Al). Then, the probability to pick 
out two Si-03 bonds and one Al-03 bond among the complete ensemble of 3xA'"o3 bonds 
can be written for the 0Si2Al triclusters, using the standard combinatory notations: 



Ai2 



Nsi-03 



^Al-03 



CI 



(7) 



3xAfo3 



where 



= C{n, k) 



k\{n-ky: 



(8) 
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The derivation of this formula for the other tricluster types is triviaL When we compare 
the expected values for a random distribution of Si and Al among the possible sites to 
the values obtained in the simulated sample, an excess of 0SiAl2 is observed in the 1600- 
atoms system since 13 0SiAl2 triclusters are found instead of 8.2, the randomly expected. 
This excess number of 0SiAl2 is accompanied by a deficit of 0Si2Al and OAI3 triclusters, 
compared to a random distribution of Si and Al. This result suggests a possible charge 
compensation role of the oxygen triclusters and could partly explain the violation of the 
Al/Al avoidance rule since Al-O-Al linkages are present in the 0SiAl2 triclusters. 

E. Rings distribution 

In silicate systems, rings are defined as loops of T-0 links (T=Si or Al) and the ring 
size distribution is a measure of the intermediate range order. In the case of aluminosilicate 
glasses, we defined a ring as the smallest loop composed of T-0 links. Fig. ^2 presents the 
ring size distributions computed for all system sizes. The distributions extend from rings of 
size 3 to rings of size 9 and they present a maximum located at 6 for all system sizes. As the 
size of the system decreases, the size of the box becomes comparable to the size of the rings 
and only few large rings can be found. However it is interesting to note that the distribution 
maxima do not shift to smaller values for the small systems, which indicates that the finite 
size effects are not significant on this quantity. This result is interesting since we know from 
the RDFs^ that correlations exist in the 1600-atoms system up to 8 A at least, which is 
much larger than half the simulation box length of the 100-atoms system. Yet the fact that 
the medium range correlations are not taken into account in the 100-atoms system does not 
seem to affect the ring size distribution. 

As the maximum of the rings distribution is 6 for pure silica^- and 5 for aluminosilicate 
(AS2)^^, we do not denote an effect of the presence of aluminum or calcium atoms in the 
rings distribution which is silica-like. 

F. Static structure factor 

In order to compare the global structure of the simulated systems with the one of real 
systems, we computed the neutron static structure factors. This was done using the relations: 
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Sn{q) 



1 



a/3 



(9) 



with 5*0/3 (g) 



N 



$^$^(exp(iq- (r/ -r„))) 



(10) 



i=l m=l 



where q is the scattering vector, the number of atoms, Na and Nf^ the number of atoms 
of species a and /?, respectively. The factor fa/3 is equal to 0.5 for a 7^ /? and equal to 1.0 
for a = P, ba and 6/3 are the neutron scattering lengths of species a and P, and are equal to 
4.149, 5.803, 3.449, 4.700 fm for Si, O, Al and Ca, respectively^. 

Figure IT^ shows the comparison of the neutron structure factors calculated for the systems 
of 100 and 1600 atoms. Even if the 100-atoms ^^(q) is noisier, there are no perceptible 
differences between the 100- and the 1600-atoms structure factors. From the partial structure 
factors, the principal contributions to the ^^(q) first peak appear to be due to the Si-Si, 
Si-0, Al-0 and 0-0 correlations (see Fig. IT!?|l . Concerning the second peak, the main 
contributions can be attributed to Si-Si, Si-Al, 0-0 and Ca-0 correlations. It is, however, 
more difficult to make a clear relation between the other peaks and the partial correlations. 

On the other hand, in the 1600-atoms 5'„(g), one can notice the existence of a very small 
pre-peak around 0.7 (denoted by a vertical dashed line in Fig. IT^ which appears to be 
due to Si-Si, Si-0 and 0-0 correlations. Since the Al-Al correlations do not contribute to 
this pre-peak, it is very likely that this feature is indeed due to the presence of the Ca atoms 
which induce the existence of pockets or channels and consequently, the oorganisation of the 
matrix on large scales. This type of long-range organisation has already been observed in 
sodo-silicate glasses and melts^. 

Since no experimental data were available for the specific composition of the present 
work, in order to compare the computed ^^(q) with experimental data, we had to gener- 
ate a glass with a slightly different composition (60% Si02, 10% AI2O3, 30% CaO). This 
new simulation was performed on a system of 800 atoms following the same procedure as 
described beforehand. The comparison between the computed and experimental structure 
factors from Ref.'Ai is presented in Fig. The agreement is good both for the positions 
of the peaks and for their intensities except at large q where irregularities do not permit to 
superpose the shoulder after the fourth peak. 

The difference between the two computed compositions is the Si02/CaO ratio, which is 
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lower in the case of the experimental comparison. The differences between the structure 
factors at the two compositions, Fig. and Fig. El mainly concern the two first peaks. 
In the experimental composition the first peak is lower and the second one higher than 
in the initial studied composition. This shift is due to the larger quantities of Ca in the 
experimental composition and has already been observed in a neutron scattering study of 
calcium aluminosilicate glasses of different compositionsi^. 

V. DISCUSSION 

The set of interatomic potentials used to simulate the calcium aluminosilicate glasses 
leads to structures that reproduce a series of experimental features if the studied sample 
contains at least 200 atoms. 

Firstly, the simulated neutron structure factor is globally in agreement with the exper- 
imental one although some of the first neighbor distances are not correctly reproduced. 
Recent reverse Monte Carlo calculations^^ carried out on glasses of several compositions 
have suggested that the distance discrepancies in the local ionic environments could ex- 
plain the slight differences between the simulated and experimental structure factors. By 
stretching the local Si-0 distances by 0.02 A and reducing the Al-0 and Ca-0 distances by 
respectively 0.05 A and 0.1 A, these differences could be removed. 

Secondly, the Al/Al avoidance rule is qualitatively reproduced but we have no quantitative 
experimental results to compare with simulations. The fact that the avoidance rule is not 
completely obeyed has already been suggested in recent NMR experiments^. The authors 
predict that there still exist Al-O-Al linkages in a sodium aluminosilicate glass with Si/Al = 
0.7 and in a calcium aluminosilicate glass with Si/Al = 0.5 which suggests that the respective 
concentrations in Si-O-Si, Si-O-Al and Al-O-Al are dependent on the energetic cost of the 
local configurations surrounding each linkage. 

Finally, another interesting structural feature concerns the existence of oxygen triclusters. 
If we consider that each Ca atom is "used" to compensate the charged AIO4 tetrahedra and 
that the remaining Ca ions create non-bridging oxygens, the predicted concentration of non- 
bridging oxygens is lower than the simulated ones. This result means that the Ca atoms 
are not completely consumed by the AIO4 tetrahedra, and part of the AIO4 tetrahedra are 
charge-balanced by the oxygen triclusters. The existence of the latter, proposed originally 
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in Ref.— , has not been directly observed but is deduced from ^ O NMR measurements 
in aluminosilicate glasses. The existence of 5 % of non-bridging oxygen in a CaAl2Si208 
glass, detected by ^'^O NMR spectroscopy, suggests also the presence of oxygen triclustera^. 
Previous simulations based on similar potentials have also predicted the existence of oxygen 
triclusters in a CaAl2Si208 glass^*^. 

However for the glass models containing 100 atoms, the structural characteristics are quite 
different indicating that these systems are subject to finite size effects. To test whether these 
size effects are due to the used potential or to the complexity of the glass, we performed 
simulations on a 78-atoms system of pure a-Si02 by using two different potentials: the 
present one detailed in section IH] and the so-called BKS potential derived by van Beest 
a/.—. Whereas the Si02 glass obtained with the BKS potential did not show any finite size 
effects, the size effects observed in the present work on the RDF first peaks and on the 
angular distributions were still present with the potential described in section [n] The main 
differences between the two potentials reside in the fact that the BKS potential comprises 
partial charges for the Si and O ions and does not include three body terms. 

In order to determine whether the size effects in the CAS system are induced by the 
presence of the three body terms, we have performed calculations without the three body 
terms. By neglecting the angular terms, we do not obtain a correct structure, since a 
large number of miscoordinated silicon atoms are found. However the structural quantities 
concerned by the finite size effects can still be examined. For instance, the first peaks of the 
Si-0 and Al-0 radial distribution functions for systems of 100 and 1600 atoms have been 
compared without the three body terms and we clearly observe the disappearance of the 
size effects on those RDFs. The finite size effects also disappear in the angular distributions. 
These considerations allow us to clarify the effect of the angular terms on the calculated 
glassy structures. The enhancement of the local order observed for the small systems is due 
to the three body terms: it completely disappears when these terms are removed. 

However, for the smallest systems, finite size effects are still present since the structure on 
distances beyond L/2 is influenced by the periodic boundary conditions. Therefore it is not 
possible to investigate the intermediate range order, which in these type of materials extends 
up to 10-15 A, in systems that are smaller than 10 A. These effects are not responsible of 
the enhancement of the local structure, but rather of the modification of the structure factor 
at low q. These effects are somehow trivial and can be easily be taken into account. 
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VI. CONCLUSION 



In this work, we have presented a systematic study of the hquid dynamics and the struc- 
tural properties of a calcium aluminosilicate glass as a function of the system size. This 
study was carried out by means of molecular dynamics simulations using a potential made 
of two-body and three-body interactions. The results show that finite size effects exist for 
the smallest system size (100 atoms) and affect a large part of structural characteristics: the 
Si-O-Si and Si-O-Al angular distributions, the Ca coordination numbers, the first peaks of 
the Si-0, Al-0 and Ca-0 pair distribution functions, and the relative number of bridging 
oxygens, non-bridging oxygens and oxygens triclusters. 

Most of these size effects were clearly attributed to the presence of three-body interactions 
in the potential. Interestingly, no finite size effects were visible in the relative percentages of 
Al-O-Al linkages and on the ring size distribution. For the larger system size, the structural 
characteristics were compared to experimental data when possible and a correct agreement 
was obtained for the Si-0 and Al-0 interatomic distances and for the neutron static structure 
factors. However, the Si-O-Si and Si-O-Al mean angles and the Ca coordination numbers 
were found to be slightly too large. 

In conclusion, the present work gives confidence in the fact that the local structure of 
relatively complex glasses can be studied by the means of ab initio simulations carried out 
on relatively small samples (200 atoms). On the other hand, it also gives indications on the 
weaknesses of the presently used classical potential and guidelines for possible improvements. 
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TABLE I: Parameters of the two-body potentials. 



atomic pair Aij [eV] 5^ [A] Qj [eV] Dij [eV] 

Si-Si 875.9321 3.448 

Si-0 1044.9519 3.049 

Si-Al 956.3850 3.448 

Si-Ca 4000.4894 3.448 

0-0 369.2783 2.857 

0-Al 1725.0881 3.448 

0-Ca 8973.4045 3.448 539.9580 884.2333 

Al-Al 1039.5842 3.448 

Al-Ca 4326.2957 3.448 

Ca-Ca 17897.5012 3.448 



TABLE II: Parameters of the three-body potentials. 



atomic triplet 


Ajifc [eV] 


7ii [A] 


lik [A] 


4 [A] 


4 [A] 


djik [degree] 


0-Si-O 


149.7960 


2.6 


2.6 


3.0 


3.0 


109.5 


0-Al-O 


149.7960 


2.6 


2.6 


3.0 


3.0 


109.5 


Si-O-Si 


6.2415 


2.0 


2.0 


2.6 


2.6 


160.0 
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TABLE III: Interatomic 


distances in the CAS glass 


for a system of 1600 atoms at 300 K, as 


extracted from the positions of the first peak maxima 
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Al-Al 
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TABLE IV: Percentages of the Si tetrahedra in the environment for all system sizes. 





100 


200 


System size 
400 


800 


1600 


Qi 





0.6 





0.6 





Q2 


2.9 


1.7 


2.3 


3.4 


3.7 


Qs 


22.7 


23.8 


26.1 


23.9 


25.3 


Q4 


74.4 


73.9 


71.6 


72.1 


71.0 
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FIGURE CAPTIONS 

Figure 1 : Time and temperature dependence of the mean squared displacements of Si, 
O, Al and Ca atoms for the 200-atoms systems at 3300 K and 5800 K. 

Figure 2 : Inverse temperature dependence of the self diffusion constants for Si, O, Al, 
and Ca atoms in the 200-atoms liquids. 

Figure 3 : (a) Si-0, Al-0 and Ca-0 radial distribution functions for systems of 100 
(dashed hues), 200, 400, 800 and 1600 atoms (bold hues), (b) Zoom on the first peak of the 
Si-0 radial distribution function, (c) Zoom on the first peak of the Al-0 radial distribution 
function, (d) Zoom on the first peak of the Ca-0 radial distribution function. 

Figure 4 : (a) Si-Ca, (b) Ca-Ca, (c) Al-Ca and (d) Al-Al radial distribution functions 
for the systems of 100 (dashed lines) and 1600 atoms (bold fines). 

Figure 5 : (a) Si-0 radial distribution functions for different Si environments Q2, Q3, Q4 
and for all Si atoms, (b) Ca-0 radial distribution functions for different oxygen atom types 
BO, NBO and for all oxygen atoms, (c) Si-Ca radial distribution functions for different Si 
environments Q2, Q3, Q4 and for all Si atoms. All RDFs are computed for the 1600-atoms 
system. 

Figure 6 : Angular distributions of the (a) Al-O-Al, (b) Si-O-Si, (c) Si-O-Al, (d) 0-Al-O 
and (e) 0-Si-O angles for systems of 100 (dashed lines) and 1600 (bold lines) atoms. 

Figure 7 : Number of linkages divided by the number of oxygen atoms as the function 
of the system size: (a) Si-O-Si, (b) Si-O-Al, (c) Al-O-Al and (d) sum of Si-O-Si, Si-O-Al 
and Al-O-Al. 

Figure 8 : Averaged probability distribution for finding N oxygen atoms around the 
calcium atoms in a system of 1600 atoms and for different cutoff radii. The inset shows 
schematically how the cutoff radii were chosen around the minimum after the first peak of 
the Ca-0 radial distribution function. 
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Figure 9 : Probability distribution for finding N oxygen atoms around the calcium atoms 
for systems from 100 to 1600 atoms. 

Figure 10 : Concentration of (a) NBO atoms, (b) BO atoms and (c) oxygen triclusters 
as the function of the system size. 

Figure 11 : Probability of finding rings of N tetrahedra for systems of 100, 200, 400, 800 
and 1600 atoms. 

Figure 12 : Static neutron structure factors for 100-atoms (bold line) and 1600-atoms 
(circles) systems. The scattering lengths were taken to be equal to bsi,o,Ai,Ca = 4.149, 5.803, 
3.449, 4.700 fm. 

Figure 13 : Partial structure factors for the 1600-atoms system. The curves have been 
shifted upward and downward for clarity and the vertical dashed line denotes the position 
of the pre-peak. 

Figure 14 : Comparison of the experimental (from Ref.J^) and calculated neutron 
structure factors for a 800-atoms CAS glass of composition 60% Si02 - 10% AI2O3 - 30% 
CaO. The scattering lengths were taken to be equal to bsi,o,Ai,Ca = 4.149, 5.803, 3.449, 
4.700 fm. 
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